/*
 *  cylk_imageManagement.c
 *  cylk
 *
 *  Created by Benjamin Heasly on 7/14/09.
 *  Copyright 2009 __MyCompanyName__. All rights reserved.
 *
 */

#include "cylk_imageManagement.h"

#pragma mark cylkImage
cylkImage new_cylkImage(int nx, int ny, int nz, float dx, float dy, float dz) {
	fprintf(stderr, "---- make image with size [%d %d %d] with scales [%.3f %.3f %.3f]\n",
			nx, ny, nz, dx, dy, dz);

    cylkImage im;
    im.nx = nx;
    im.ny = ny;
    im.nz = nz;
    im.nVoxels = nx*ny*nz;
    im.dx = dx;
    im.dy = dy;
    im.dz = dz;
    im.pt = calloc(im.nVoxels, sizeof(*(im.pt)));
    return(im);
}

cylkImage new_cylkImageLikeExisting(cylkImage *existing) {
    cylkImage im = new_cylkImage(existing->nx, 
								 existing->ny, 
								 existing->nz,
								 existing->dx, 
								 existing->dy, 
								 existing->dz);
    return(im);
}

void destroy_cylkImage(cylkImage *im) {
	free(im->pt);
	im->pt = NULL;
}


#pragma mark cylkSlice
cylkSlice new_cylkSlice(int nx, int ny, float dx, float dy) {
    cylkSlice sl;
    sl.nx = nx;
    sl.ny = ny;
    sl.nVoxels = nx*ny;
    sl.dx = dx;
    sl.dy = dy;
    sl.pt = calloc(sl.nVoxels, sizeof(*(sl.pt)));
    return(sl);
}

cylkSlice new_cylkSliceFromCylkImage(cylkImage *im, int slicenum, int slicedim) {
	cylkSlice sl;
	long start_ind, block_size, cursor_loc;
	int j, k;
	switch (slicedim) {
		case 2:
			/* extract a z slice */
			sl = new_cylkSlice(im->nx, im->ny, im->dx, im->dy);
			block_size = sl.nx * sl.ny;
			subscriptsToIndex(im, 0, 0, slicenum, &start_ind);
			memcpy(sl.pt, (im->pt + start_ind), block_size*sizeof(short));
			break;
		case 1:
			/* extract a y slice */
			sl = new_cylkSlice(im->nx, im->nz, im->dx, im->dz);
			cursor_loc = 0;
			block_size = sl.nx;
			for (k=0; k<im->nz; k++) {
				subscriptsToIndex(im, 0, slicenum, k, &start_ind);
				memcpy((sl.pt + cursor_loc), (im->pt + start_ind), block_size*sizeof(short));
				cursor_loc += block_size;
			}
			break;
		case 0:
			/* extract an x slice */
			sl = new_cylkSlice(im->ny, im->nz, im->dy, im->dz);
			for (k=0; k<im->nz; k++) {
				for (j=0; j<im->ny; j++) {
					subscriptsToIndex(im, slicenum, j, k, &start_ind);
					sl.pt[k*sl.nx + j] = im->pt[start_ind];
				}
			}
			break;
	}
	return (sl);
}

void destroy_cylkSlice(cylkSlice *sl) {
	free(sl->pt);
	sl->pt = NULL;
}


#pragma mark indexing
void subscriptsToIndex(cylkImage *im, int x, int y, int z, long *index) {
	*index = x + y*im->nx + z*im->nx*im->ny;
}

void indexToSubscripts(cylkImage *im, int *x, int *y, int *z, long index) {
	*x = index % im->nx;
    *y = (index / im->nx) % im->ny;
    *z = index / (im->nx*im->ny);
}

#pragma mark misc

long indexOfMaxImageVoxel(cylkImage *im) {
	long i, maxIndex=0;
	for (i=0; i<im->nVoxels; i++) {
		if (im->pt[i] > im->pt[maxIndex])
			maxIndex = i;
	}
	return(maxIndex);
}

int indexesOfVoxelsNeighboringImageVoxel(cylkImage *im, long voxelIndex, long *neighborIndexes) {
    int i, x, y, z, xNeighbor, yNeighbor, zNeighbor, numNeighbors=0;
    long neighborIndex;
	
	// all possible neighbor voxels pixels, including diagonals
	int neighbors[26][3] = {
		{-1,-1,-1},  {-1,-1,0}, {-1,-1,1}, {-1,0,-1}, {-1,0,0}, {-1,0,1}, {-1,1,-1}, {-1,1,0}, {-1,1,1},
		{ 0,-1,-1},  { 0,-1,0}, { 0,-1,1}, { 0,0,-1},           { 0,0,1}, { 0,1,-1}, { 0,1,0}, { 0,1,1},
		{ 1,-1,-1},  { 1,-1,0}, { 1,-1,1}, { 1,0,-1}, { 1,0,0}, { 1,0,1}, { 1,1,-1}, { 1,1,0}, { 1,1,1}};
	
	// get subscripts form voxel index
    indexToSubscripts(im, &x, &y, &z, voxelIndex);
    
    for(i=0; i<26; i++) {
        xNeighbor = x + neighbors[i][0];
        yNeighbor = y + neighbors[i][1];
        zNeighbor = z + neighbors[i][2];
		
		// add all neighboring pixels that are in image bounds 
        if ((xNeighbor>=0) && (xNeighbor<im->nx) && (yNeighbor>=0) && (yNeighbor<im->ny) && (zNeighbor>=0) && (zNeighbor<im->nz)) {
			subscriptsToIndex(im, xNeighbor, yNeighbor, zNeighbor, &neighborIndex);
            neighborIndexes[numNeighbors++] = neighborIndex;
        }
    }
	return(numNeighbors);
}

int insertSliceIntoImage(cylkSlice *inputSlice, cylkImage *outputImage, int slicenum,int slicedim) {
	long start_ind, cursor_loc, block_size;
	int j, k;
	switch (slicedim) {
		case 2:
			//insert a z slice
			if (inputSlice->nx != outputImage->nx
				|| inputSlice->ny != outputImage->ny)
				return (1);
			block_size = inputSlice->nx * inputSlice->ny;
			subscriptsToIndex(outputImage, 0, 0, slicenum, &start_ind);
			memcpy ((outputImage->pt + start_ind), inputSlice->pt,block_size * sizeof (short));
			break;
		case 1:
			// insert a y slice
			if (inputSlice->nx != outputImage->nx || inputSlice->ny != outputImage->nz)
				return (1);
			cursor_loc = 0;
			block_size = inputSlice->nx;
			for (k = 0; k < inputSlice->ny; k++) {
				subscriptsToIndex(outputImage, 0, slicenum, k, &start_ind);
				memcpy ((outputImage->pt + start_ind),(inputSlice->pt + cursor_loc),sizeof (short) * block_size);
				cursor_loc += block_size;
			}
			break;
		case 0:
			// insert an x slice
			if (inputSlice->nx != outputImage->ny
				|| inputSlice->ny != outputImage->nz)
				return (1);
			for (k = 0; k < inputSlice->ny; k++) {
				for (j = 0; j < inputSlice->nx; j++) {
					subscriptsToIndex(outputImage, slicenum, j, k, &start_ind);
					outputImage->pt[start_ind] =
					inputSlice->pt[k * inputSlice->nx + j];
				}
			}
			break;
	}
	return (0);
}